#include <fintrf.h>
!
! API here: http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/apiref/bqoqnz0.html&http://www.mathworks.com/support/tech-notes/1600/1605.html
!
      SUBROUTINE mexFunction(NLHS, PLHS, NRHS, PRHS)
      implicit none;

      ! Pointers for bringing in in/out arguments
      mwPointer PLHS(*), PRHS(*)
      mwsize NLHS, NRHS

      ! Define the functions to use
      mwPointer mxCreateDoubleMatrix, mxGetPr

      ! Functions to get matrix sizes
      mwSize mxGetM, mxGetN

      ! Testing functions
      INTEGER*4 mxIsNumeric

      ! Used to size matrixes
      mwSize m, n, size

      ! Output pointers
      mwPointer x_p
      mwPointer y_p, yd_p
      mwPointer u_p, ud_p
      mwPointer dx_p, dt_p

      ! A character array to use for writing output messages to matlab
      character*120 :: line

      ! Variable to temporarily store input variables as reals before they are converted to ints
      REAL*8 :: tmp_real

      ! The actual program variables

      ! Define all my parameters
      integer*4 :: dx = 1000 ! [m]
      integer*4 :: dt = 240; ! [s]
      INTEGER*4 :: max_signal_length, i;
      parameter (max_signal_length=100)

      ! Data
      DOUBLE PRECISION, DIMENSION(max_signal_length) :: x0, x, y
      DOUBLE PRECISION, DIMENSION(max_signal_length) :: xd

      ! Wind
      DOUBLE PRECISION, DIMENSION(max_signal_length) :: u
      DOUBLE PRECISION, DIMENSION(max_signal_length) :: ud

      ! Diffusion
      DOUBLE PRECISION, DIMENSION(max_signal_length) :: k
      DOUBLE PRECISION, DIMENSION(max_signal_length) :: kd

      !
      ! Initialize the data
      !

      ! Wind
      data u /max_signal_length*1.0/ ! [m/s]

      ! Diffusion
      data k /max_signal_length*100/ ! [m^2/s]

      ! Initialize derivatives (SHOULDN'T HAVE TO DO THIS.. Hopefully
      ! I can remove this soon
      data ud /max_signal_length*0.0/
      data kd /max_signal_length*0.0/

      ! CHECK FOR PROPER NUMBER OF ARGUMENTS

      IF (NRHS .NE. 4) THEN
        CALL mexErrMsgTxt('advf requires four inputs: [x u dx dt]')
      ELSEIF (NLHS .NE. 2) THEN
        CALL mexErrMSGTxt('Outputs: [y yd]')
      ENDIF

      ! Get the size of the input array
      m = mxGetM(prhs(1))
      n = mxGetN(prhs(1))
      size = m*n

      ! Make sure the input isn't exceeding constraints
      if (m .gt. 1) then
        call mexErrMsgTxt("X input must me a vector")
      elseif (size .gt. max_signal_length) then
        call mexErrMsgTxt("X vector must be <= max_signal_length elements in length")
      endif

      ! Ensure the array is numeric (not strings)
      if (mxIsNumeric(prhs(1)) .eq. 0) then
        call mexErrMsgTxt('X vector must be numeric')
      endif

      ! Create an output matrix for the AD'ed signal
      PLHS(1) = mxCreateDoubleMatrix(m,n,size)
      PLHS(2) = mxCreateDoubleMatrix(m,n,size)

      ! Assign pointers to the various parameters
      ! Output
      y_p      = mxGetPr(PLHS(1))
      yd_p     = mxGetPr(PLHS(2))

      ! Input
      x_p      = mxGetPr(PRHS(1))
      CALL mxCopyPtrToReal8(x_p, x, size)

      u_p      = mxGetPr(PRHS(2))
      CALL mxCopyPtrToReal8(u_p, u, size)

      ud_p     = mxGetPr(PRHS(3))
      CALL mxCopyPtrToReal8(ud_p, ud, size)

      dx_p      = mxGetPr(PRHS(4))
      CALL mxCopyPtrToReal8(dx_p, tmp_real, int(1, 4))
      dx = int(tmp_real, 4)

      dt_p      = mxGetPr(PRHS(5))
      CALL mxCopyPtrToReal8(dt_p, tmp_real, int(1, 4))
      dt = int(tmp_real, 4)

      !
      ! Copy x to x0
      !
      x0 = x;

      !
      ! Perform the actual computations..
      !

      CALL advf_d(x0, x, xd, u, k, int(size, 4), dx, dt)

      !write(line,*) x(7:13), "\n"
      !CALL mexPrintf(line)

      ! COPY OUTPUT WHICH IS STORED IN LOCAL ARRAY TO MATRIX OUTPUT
      CALL mxCopyReal8ToPtr(x, y_p, size)
      CALL mxCopyReal8ToPtr(xd, yd_p, size)

      RETURN
      END

      ! Include our subroutines
      include "advf_d.f95"

